data{
  int<lower=0> N;
  int<lower=1> K; //number of category
  int<lower=0, upper=1> treatment[N];
  int<lower=1,upper=K> Y[N];
}
parameters{
  real beta; // beta for treatment
  ordered[K-1] c;
}
model{
  for(i in 1:N){
    Y[i] ~ ordered_logistic(beta*treatment[i], c);
  }

}
generated quantities{
	real T[K];
	real C[K];
	real diff[K];
	real temp1;
	real temp2;

  real diff_D[2];  // dichotomized
  real diff_DwN[2];  // dichotomoized with neutral

	for(k in 1:K){
		T[k]=0;
		C[k]=0;
		diff[k]=0;
	}

  for(k in 1:2) {
    diff_D[k] = 0;
    diff_DwN[k] = 0;
  }

	for(i in 1:N){

		// See a node in ordered_logit3.stan

		T[1] = inv_logit(- beta*1 + c[1]);
		C[1] = inv_logit(- beta*0 + c[1]);	
		diff[1] = diff[1] + T[1] - C[1];

		for(s in 2:K-1){

			temp1 = inv_logit(- beta*1 + c[s-1]);
			temp2 = inv_logit(- beta*1 + c[s]);
			T[s] =  fmax(temp1, temp2) - fmin(temp1, temp2);


			temp1 = inv_logit(- beta*0 + c[s-1]);
			temp2 = inv_logit(- beta*0 + c[s]);
			C[s] =  fmax(temp1, temp2) - fmin(temp1, temp2);
			diff[s] = diff[s] + T[s] - C[s];
		}

		T[K] =  1 - inv_logit(- beta*1 + c[K-1]);
		C[K] =  1 - inv_logit(- beta*0 + c[K-1]);
		diff[K] = diff[K] + T[K] - C[K];

    diff_D[1] = diff[1] + diff[2];
    diff_D[2] = diff[4] + diff[5];

    diff_DwN[1] = diff[1] + diff[2];
    diff_DwN[2] = diff[3] + diff[4] + diff[5];

	}

	for(k in 1:K) {
		diff[k] = diff[k]/N;
  }

  for (k in 1:2) {
    diff_D[k] = diff_D[k]/N; 
    diff_DwN[k] = diff_DwN[k]/N; 
  }

}

